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Abstract: A symmetric discretisation scheme for heterogeneous anisotropic diffusion problems on general 
meshes is developed and studied. The unknowns of this scheme are the values at the centre of the control 
OO ' volumes and at some internal interfaces which may for instance be chosen at the diffusion tensor discontinuities. 

The scheme is therefore completely cell-centred if no edge unknown is kept. It is shown to be accurate on several 
numerical examples. Convergence of the approximate solution to the continuous solution is proved for general 
(possibly discontinuous) tensors, general (possibly nonconforming) meshes, and with no regularity assumption 
on the solution. An error estimate is then deduced under suitable regularity assumptions on the solution. 

' Keywords : Heterogeneous anisotropic diffusion, nonconforming grids, finite volume schemes 

1 Introduction 

< 

, Anisotropic heterogeneous diffusion problems arise in a wide range of scientific fields such as hy- 

drogeology, oil reservoir simulation, plasma physics, semiconductor modelling, biology, etc.. When 
implementing numerical methods for this kind of problem, one needs to find an approximation of u, 
weak solution to the following equation: 



-div(A(a;)V'u) = / in 17, (1) 

in 

^ i with boundary condition 

ti = on dn, (2) 



O 



where we denote by dQ = Q\Q the boundary of the domain 0, under the following assumptions: 

^ . is an open bounded connected polyhedral subset of M'^, d G N \ {0}, (3) 

00 . 

, A is a measurable function from 17 to A^^(]R), (4) 

where we denote by A4dO^) the set oi dx d matrices, such that for a.e. x £ ^l, ^i^) is symmetric, and 
^ . such that the set of its eigenvalues is included in [A, A], with A and A G M satisfying < A < A, and 

f€L\n). (5) 



Under these hypotheses, the weak solution of (l)-(2) is the unique function u satisfying: 

u G H^{n), 

(6) 



A{x)Vuix)-Vv{x)dx= [ f{x)v{x)dx Vt; G i/o(0). 

Ju 



Usual discretisation schemes for Problem (6) include finite difference, finite element or finite volume 
methods. Finite volume methods are actually very popular in oilreservoir engineering, a probable 
reason being that complex coupled physical phenomena may be discretised on the same grids. The 
well-known five-point scheme on rectangles (see e.g. [29]) and four-point scheme on triangles [23] are 
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not easily adapted to heterogeneous anisotropic diffusion operators [24]. A scheme with an enlarged 
stencil, which handles anisotropy on meshes satisfying an orthogonality property, was proposed and 
analysed in [18]. Another problem that has to be faced in several fields of applications (such as 
hydrogeology and oil reservoir engineering) is the fact that the discretisation meshes are imposed 
by engineering and computing considerations; therefore, we have to deal with distorted and possibly 
nonconforming meshes. 

A huge literature exists in the engineering setting, so we shall not try to be exhaustive. Let us 
nevertheless mention the finite volume schemes using the well-known multipoint flux approximation 
[1, 2, 3]. These schemes involve the reconstruction of the gradient in order to evaluate the fluxes, which 
is also the case in [13, 28]. Among other approaches let us cite [22], which uses a parametrisation 
technique. However, even though these schemes perform well in a number of cases, their convergence 
analysis often seems to remain out of reach, except under additional geometrical conditions [13]. 
More recently, flnite volume schemes using interface values have been studied. In [19] we presented 
a "hybrid finite volume" (HVF) scheme for any space dimension, which involves edge unknowns in 
addition to the usual cell unknowns, and in [15], a "mixed finite volume" scheme (MFV) was proposed, 
which involves the fluxes and the values as unknowns. This is also the case for the mimetic finite 
difference (MFD) schemes [9, 10], which were introduced previously; in spite of their name, mimetic 
schemes are very much in the finite volume spirit, since they rely on both a fiux balance equation and 
on the local conservativity of the numerical fiuxes, that are probably the two "pillars" of the finite 
volume philosophy; but then, finite volume schemes are also often called finite difference schemes in 
the engineering literature because of the finite difference approximation of the fluxes. In fact, a recent 
benchmark [25] provided sufficient information to suspect that the methods HFV, MFV and MFD 
indeed coincide at the algebraic level and establishing this is the aim of ongoing work [16]. Let us 
mention that the Raviart-Thomas mixed flnite element method, which also involves edge unknowns, 
was generalised to handle distorted hexahedral meshes [27]. These schemes require the fluxes or edge 
unknowns as additional values (or as sole values after hybridisation), and they may be more expensive 
than cell-centred schemes, especially in the 3D case. 

In the two-dimensional case, we also mention [6], which discusses a scheme based on vertex reconstruc- 
tions, and the family of double mesh schemes [26, 14, 7]. The generalisation of this type of scheme to 
3D is the subject of ongoing work. 

The scheme that we present here is designed on very general polygonal, possibly non-convex and 
nonconforming meshes, with the following two priorities in mind: 

• For cost reasons and data structure issues, we wish to obtain a symmetric scheme which is as 
close as possible to a cell-centred scheme, that is to a scheme involving one unknown per control 
volume (or grid cell). 

• For accuracy reasons, we require the local conservativity of the numerical flux to hold at the 
interfaces between highly heterogeneous media. 

In [21], we introduced a cell-centred scheme for the approximation of the Laplace operator on noncon- 
forming grids in the framework of the incompressible Navier-Stokes equations [21] and which may be 
viewed as a low order nonconforming Galerkin approximation. The scheme (called "SUCCES" in [4]) 
was also implemented for anisotropic and heterogeneous problems on general meshes, and was shown 
to be highly competitive for oil reservoir simulation in comparison with other well-known schemes 
such as the multiple point flux approximation schemes. It is cheaper than the above mentioned hybrid 
type schemes (HVF, MFV and MFD) because it is based on cell unknowns only. However, it is not 
as accurate as the hybrid schemes for strongly heterogeneous problems, very likely because of the 
weaker approximation of the normal fluxes at the heterogeneous interfaces. In the present work, we 
construct a discretisation scheme (SUSHI) for any kind of polyhedral mesh, which incorporates the 
best properties of the cell-centred (SUCCES) and hybrid (HFV) schemes: unknowns on the edges are 
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only introduced when needed, for instance when there is strong medium heterogeneity at these edges. 
If the set of edge unknowns is empty, then SUSHI reduces to the above mentioned cell-centred scheme; 
if unknowns are associated to all internal edges, then SUSHI is the hybrid scheme HFV. 

The outline of this paper is as follows. In Section 2, we present the guidelines which led us in the 
construction of convergent schemes on general nonconforming meshes. The practical properties of 
the resulting schemes are shown through numerical examples in Section 3. Then the mathematical 
analysis of convergence and error estimation are performed in Section 4. This analysis is based on 
some discrete functional analytic tools, such as discrete Sobolev inequalities, which are provided in 
Section 5. Conclusions and perspectives are discussed in Section 6. 

2 Fundamentals for a class of nonconforming schemes 

Let us first present the desired properties which have led us to the design of the schemes under study: 

(PI) The schemes must apply on any type of grid: conforming or nonconforming, 2D and 3D (or more, 
see for instance the frameworks of kinetic formulations or financial mathematics), consisting of 
control volumes which are only assumed to be polyhedral (the boundary of each control volume 
is a finite union of subsets of hyperplanes). 

(P2) The matrices of the linear systems generated are expected to be sparse, symmetric and positive 
definite. 

(P3) We wish to be able to prove the convergence of the family of discrete solutions to the solution 
of the continuous problem as the mesh size tends to 0, and of the family of associate gradients 
to the gradient of the solution, with no regularity assumption on the solution of the continuous 
problem, and to derive error estimates when the analytic solution is regular enough. 

In order to describe the schemes we now introduce some notations for the space discretisation. 

Definition 2.1 (Space discretisation) Let Q he a polyhedral open bounded connected subset ofW^, 
with d G N \ {0}, and dO, = its boundary. A discretisation ofQ, denoted by D, is defined as the 

triplet V = {M,8,V), where: 

1. M is a finite family of nonempty connected open disjoint subsets of (the "control volumes") 
such that Q. = UkgmK- For any K G M., let dK = K\K be the boundary of K; let \K\ > 
denote the measure of K and let hx denote the diameter of K . 

2. £ is a finite family of disjoint subsets of (the "edges" of the mesh), such that, for all a ^ £, 
a is a nonempty open subset of a hyperplane ofW^, whose {d — 1)- dimensional measure \a\ is 
strictly positive. We also assume that, for all K £ Ai, there exists a subset £k of £ such that 
dK = [Jfj^Sj^a. For any a £ £ , we denote by = {K £ Ai,a £ £k}- We then assume that, 
for all a £ £, either has exactly one element and then a C dVL (the set of these interfaces, 
called boundary interfaces, is denoted by £c:^t) or M.fy has exactly two elements (the set of these 
interfaces, called interior interfaces, is denoted by £mt)- For all a £ £ , we denote by Xfy the 
barycentre of a. For all K £ M and a £ £k, we denote by nK,a the unit vector normal to a 
outward to K . 

3. V is a family of points ofVL indexed by M, denoted by V = {xK)KeM; such that for all K £ M, 
xk £ K and K is assumed to be xk -star-shaped, which means that for all x £ K, the inclusion 
[xk,x\ C K holds. Denoting by dK,a the Euclidean distance between xk and the hyperplane 
including a, one assumes that dK,a > 0. We then denote by DK,a the cone with vertex xk and 
basis cr. 
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Remark 2.1 The above definition applies to a large variety of meshes. Note that no hypothesis is 
made on the convexity of the control volumes; in fact, generalised hexahedra, i.e. with faces which may 
he composed of several planar suh- faces may he used. Often encountered in subsurface flow simulations, 
such hexahedra may have up to 12 faces (resp. 24 faces) if each non planar face is composed of two 
triangles (resp. four triangles), hut only 6 neighbouring control volumes. 



2.1 From a "hybrid" finite volume scheme. . . 

The idea of the "hybrid" schemes (among them one may include the mixed finite elements, the mixed 
finite volume or the mimetic finite difference schemes) is to find an approximation to the solution of 
(l)-(2) by setting up a system of discrete equations for a family of values {{uK)KeMi (^o-)CTe^) in the 
control volumes and on the interfaces. The number of unknowns is therefore card(A^) + card(iS). 
Following the idea of the finite volume framework, Equation (1) is integrated over each control volume 
K € M, which formally gives (assuming sufficient regularity on u and A) the following balance 
equation on the control volume K: 

V ( - / A(a;)Vn(a;) • ni^,^d7(a;) ) = / f{x)dx. 

The fiux — k{x)Vu{x) ■ nK,(j^"y{x) is approximated by a function Fk,<t{u) of the values {{uk)k&m^ 
{ua-)(T&s) at the "centres" and at the interfaces of the control volumes (in all practical cases, Fk,(t{u) 
only depends on uk and all (itcr')o-'e^A-)- ^ discrete equation corresponding to (1) is then: 

V FkA^) = [ f{x)dx yK G M. (7) 

The values on the interfaces are then introduced so as to allow for a consistent approximation of the 
normal fiuxes in the case of an anisotropic operator and a general, possibly nonconforming mesh. We 
thus have card(<5) supplementary unknowns, and need card(<S) equations to ensure that the problem 
is well posed. For the boundary faces or edges, these equations are obtained by writing the discrete 
counterpart of the boundary condition (2): 

= V(j G £cxt- (8) 

Following the finite volume ideas, we may write the continuity of the discrete flux for all interior edges, 
that is to say: 

FkAu) + FL,a{u) = 0, for a G Sint such that Ma = {K, L}. (9) 
We now have card(A^) + card(i5int) unknowns and equations. 

Remark 2.2 In the case A{x) = A(a;)Id, on meshes satisfying an orthogonality condition as men- 
tioned in the introduction of this paper (this condition states the orthogonality between the line join- 
ing the centres of two neighbouring control volumes with their common interface, see [17, Defini- 
tion 9.1 p. 762]), a consistent numerical flux is obtained using the two-point formula F^^aiu) = 
^K\c\iuK — U(j)/dK,a! whcrc Xk is the average value for A in K. Then, writing (9) for all a G Sint 
such that Ma = {K, L}, we obtain Ua as a linear combination of uk o,nd ul- Plugging this expression 
into (7), we get a scheme with card{M) equations and card{M) unknowns (see [17, Section 11.1 pp. 
815-820] for more details). In the case of a rectangular (resp. triangular) mesh, this is the well-known 
five points (resp. four points) scheme with harmonic averages of the diffusion. 

With a proper choice of the expression Fx^aiu), which we shall introduce below, this scheme, first 
introduced in [19], is quite efficient for the simulation of fiuid fiow in heterogeneous media (where 
harmonic averages for A are preferred to arithmetic averages [5]) and may be shown to converge. This 
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scheme does have one drawback: since the number of unknowns is the sum of the number of control 
volumes and of interior interfaces, the resulting scheme is quite expensive (although it is sometimes 
possible to algebraically eliminate the values at the control volumes, as in the mixed hybrid finite 
element method, see [8, pp. 178-181]). 

Remark 2.3 Note that in the case of regular conforming simplices (triangles in 2D, tetrahedra in 
3D), there is an algebraic possibility to express the unknowns {ua)ae£ <is local affine combinations of 
the values {uk)k<^m ^"i^d therefore to eliminate them [31]. The idea is to remark that the linear system 
constituted by the equations (7) for all K G Ms, where Ms is the set of all simplices sharing the same 
interior vertex S, and (9) for all the interior edges such that Ai^^ C M.s, presents as many equations 
as unknowns u„, for a G Uk^Ms^k ■ Indeed, the number of edges in Uk^Ms^k such that M.a- ^ Ms 
is equal to the number of control volumes in Ms- Unfortunately, there is at this time no general result 
on the invertibility or the symmetry of the matrix of this system, and this method does not apply to 
other types of meshes than simplicial meshes. 

In order to reduce the computational cost of the scheme, we developed in [21] an idea which is in fact 
close to the finite element philosophy since we express the finite volume scheme in a weak form; to 
this end, let us first define the sets Xx> and Xx>fl where the discrete unknowns lie, that is to say: 

X'D = {V = {{vK)K&M,{Va)a&e):VK G M, G M}, (10) 

^x>,0 = G Xx> such that f o- = Vcr G <Scxt}- (11) 

Multiplying, for any v G ^d.Oj Equation (7) by the value vk of v on the control volume K and 
summing over K ^ M leads to: 



Using (9), we get the following discrete weak formulation: 



Find u G X-pfi such that: 

{u,v)f = Vk / f{x)dx, for all w G Xx),( 



(12) 



with 

{u,v)f= ^ FK,a{u){vK - Vg)- (13) 

Note that choosing v G Xjyfi such that = 1, = for any L G M,L ^ K and Uo- = for 
any a & £ yields (7). Similarly, choosing v G Xjyfi such that vk = ^ for any K ^ Ai, and = I 
and = for any t ^ £,t ^ a leads to (9). Therefore the hybrid finite volume scheme (7)-^(9) is 
equivalent to the discrete weak formulation (12). 



2.2 ... to a nonconforming finite element scheme. . . 

We may then choose to use the weak discrete form (13) as an approximation of the bilinear form 
a(-, •), but with a space of dimension smaller than that of Xd^q. This can be achieved by expressing 
the value of u on any interior interface a G <Sint as a consistent barycentric combination of the values 

UK'- 
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where (P^) ksm is a family of real numbers, with 7^ only for some control volumes K close to 
a, and such that 

5^ /3f = 1 and a;, = ^ /?f a;,^ Va G ^Ti^t. (15) 

This ensures that if (/? is a regular function, then = Z^i^eA^ (^afi^K) is a consistent approximation 
of ip{X(^) for cr G Si^f We recall that the values u^,a £ Sext are set to in order to respect the 
boundary conditions (2). Hence the new scheme reads: 

Find u £ Xx>fi such that u^- = Pa'^K Vo" G Si^t, and 

^ /■ ""'^ ^ (16) 

{u,v)f = /J vk / /(aj)da;, for all f G X^j^o with = /J P^vk Vcr G <Sint- 

This method has been shown in [21] to be efficient in the case of a problem where A = Id (for the 
approximation of the viscous terms in the Navier-Stokes problem). With an appropriate choice for the 
expression of the numerical flux, it also yields conservativity in a certain sense (more on this below), 
but no longer to the classical (in the finite volume framework) equation (9): indeed, since the degrees 
of freedom on the edges are no longer present, one may not use = 1 to recover (9). Note also 
that taking vk = 1 does not yield (7). This scheme has been implemented for the discretisation of 
the diffusive term in the incompressible Navier Stokes equations on general two- or three-dimensional 
grids, and gives excellent results [11, 12]. Unfortunately, because of poor approximation of the local 
flux at strongly heterogeneous interfaces, this approach is not sufficient to provide accurate results 
for some types of flows in heterogeneous media, as we shall show in Section 3. This is especially true 
when using coarse meshes, as is often the case in industrial problems. 



2.3 ... to an optimal compromise? 

Therefore we now propose a scheme which has the advantage of both techniques: we shall use equation 
(13) and keep the unknowns Ua- on the edges which require them, for instance those where the matrix A 
is discontinuous: hence (9) will hold for all edges associated to these unknowns; for all other interfaces, 
we shall impose the values of u using (14), and therefore eliminate these unknowns. Let us decompose 
the set £\at of interfaces into two nonintersecting subsets, that is: E^t = B [J7i,TC = fint \ S. The 
interface unknowns associated with B will be computed by using the barycentric formula (14). 

Remark 2.4 Note that, although the accuracy of the scheme is increased in practice when the points 
where the matrix A is discontinuous are located within the set Uo-eH ^> such a property is not needed 
in the mathematical study of the scheme. 

Let us introduce the space Xx>,b C Xx>fi defined by: 

Xv,l3 = {v € Xx> such that v^^ = for all a G Scxt and Va- satisfying (14) for all a G B}. (17) 

The composite scheme which we consider in this work reads: 

J Find u G Xxi^b such that: , . 

\ {u,v)f = T.KeM^K Jk fix)^^^ for all t; G Xx),B. 

We therefore obtain a symmetric scheme with card(Al) + card('H) equations and unknowns. It is 
thus less expensive while it remains accurate (for the choice of numerical flux given below) even in the 
case of strong heterogeneity (see section 3). 
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Note that with the present scheme, (9) holds for ah a € TC, but not generally for any a € B. However, 
fluxes between pahs of control volumes can nevertheless be identified. Indeed, we may write 



and therefore: 



where 



and 



{u,v)f= Y ^ FK,aiu){vK -Va) + ^ ^ Fk,l{u){vk - Vl) , 
K&M o-G^xnW {-ft',L)eA^-D 



Mv = {{K, L) G M^ 3a G 8k n i3, ^ or 3a G £:l n B, (3^ / 0} 



dK 



crefxriB o-G^inB 

Note that, if {K,L) G Mv, then {L,K) G Nv and Fk,l{u) = -Fl,k{u); furthermore, Fk,l{u) / 
implies {K,L) G Mv, and the scheme's stencil is determined by the set {L E M. such that {K,L) G 
Mv}- Then, taking vk = ^ and all other degrees of freedom of v G Xv,b equal to 0, (18) yields 



{K,L)&J\rr, 

which shows the "finite volume philosophy" of the scheme. 

Remark 2.5 (Other boundary conditions) In the case of Neumann or Robin boundary condi- 
tions, the discrete space Xv,b is modified to include the unknowns associated to the corresponding 
edges, and the resulting discrete weak formulation is then straightforward. 

Remark 2.6 (Extension of the scheme) For consistency reasons, it is preferable that the coeffi- 
cients j3^ associated with a ^ B be nonzero for points xk that lie in the same "regularity zone" of 
the solution as x„ (that is with a zone with no diffusion tensor discontinuity). This is not always 
easy: indeed, in the tilted barrier example described in Section 3.3 below, the barrier contains only one 
layer of grid cells, so that, for an internal interface of this layer, it is difficult to use points xl that 
are located in the same diffusion regularity zone with respect to xk- There is, however, no additional 
difficulty to replace (14) in the definition of (17) by 

U^= J2 P^^K + Y ^ ^' (19) 

/3f + J] /3^' = 1 and = P^^K + Z^-'^-' ^ (20) 

K&M a'€H KeM cr'en 

This trick solves the consistency issue without switching the edge to the hybrid set 7i, while all the 
mathematical properties shown below still hold. 
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2.4 Construction of the fluxes using a discrete gradient 



For the definition of the schemes to be complete, there now remains to explain how we find a convenient 
expression for Fk,(t{u) with respect to the discrete unknowns. An idea that has been used in several 
of the schemes referred to in the Introduction is to look for a consistent expression of the flux by using 
adequate linear combinations of the unknowns; however, referring to the beginning of Section 2, such 
a reconstruction does not in general lead to the desired properties (P2) (symmetric definite positive 
matrices) and (P3) (convergence). Our idea here is different: it is based on the identification of the 
numerical fluxes FK,a{u) through the mesh-dependent bilinear form {■,-)f defined in (13), using the 
expression of a discrete gradient. Indeed let us assume that, for all u £ Xx>, we have constructed a 
discrete gradient Vx)U, we then seek a family {FK,a{u)) ^eM such that 

{u,v)f= ^ ^ FK,a{u){vK - Va) = / Vt)u{x) ■ A{x)\7T>v{x)dx \/u,v G Xj,. (21) 

Remark 2.7 (On the construction of the discrete fluxes) Note that it is always possible to de- 
duce an expression for FK,a{u) satisfying (21), under the sufficient condition that, for all K ^ M and 
a.e. X £ K , Vxiu{x) is expressed as a linear combination of {u^ — UK)a^£K! coefficients of which 
are measurable bounded functions of x. This property is ensured in the construction ofVx>u{x) given 
below. 

We prove in Section 4 below that the desired properties (P2) and (P3) hold if the discrete gradient 
satisfies the following properties: 

1. (Weak compactness) For a sequence of space discretisations of O with mesh size tending to 0, if 
the sequence of associated grid functions is bounded in some sense, then their discrete gradient 
converges at least weakly in L'^iyiY to the gradient of an element of Hq{Q); 

2. (Consistency) If is a regular function from O to M, the discrete gradient of the piece-wise 
function defined by taking the value ip{xK) on each control volume K and '^{X(j) on each edge 
(T is a consistent approximation of the gradient of ip. 

Let us first define: 

^KU = \a\{ua - UK)nK,cj \/K e M,\fu e Xv, (22) 
\K\ ^ — ^ 

where nx,o- is the outward to K normal unit vector, and |cr| are the usual measures (volumes, 
areas, or lengths) of K and a. The consistency of formula (22) stems from the following geometrical 
relation: 

^\a\nK,a{xa-XK)^ = \K\\d MK^M, (23) 

where {Xfj — x^Y is the transpose of x^^ — xk G IR'^, and Id is the d x d identity matrix. Indeed, 
for any linear function defined on Q by il^{x) = G ■ x with G G M*^, assuming that Ug- = "0(3^0-) and 
UK = V'(*/<)) we get U(j — UK = {x^- — xkYG = {x^- — xkYV^, hence (22) leads to Vku = S/ip. 
Since the coefficient of uk in (22) is in fact equal to zero, a re-construction of the discrete gradient 
Vdu solely based on (22) cannot lead to a definite discrete bilinear form in the general case. Hence, 
we now introduce a stabilised gradient: 

V/<,<jM = Vku + RK,aU nK,a, (24) 

with 

RK,aU=- {u„ - UK -Vku ■ {x^ - xk)) , (25) 

(J'K^a 
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(recall that d is the space dimension and dK,a is the Euclidean distance between and a). We may 
then define Vxih as the piece-wise constant function equal to V K,aU a.e. in the cone D^^a with vertex 
Xk and basis a: 

Vt>u{x) = V K,aU for a.e. x € Dk,^- (26) 
Note that, from the definition (25), thanks to (23) and to the definition (22), we get that 

Y.^-^^^RK,aunK,a = ^ yKeM. (27) 

We prove in Lemmata 4.2 and 4.3 below that the discrete gradient defined by (22)-(26) indeed satisfies 
the above stated weak compactness and consistency properties. In order to identify the numerical 
fluxes Fk,(t{u) through Relation (21), we put the discrete gradient in the form 



fjfj 



with 



\a\ 


nK,a + 


Vd 1 


\K 


dK,cr \ 


W\ 






\K 




dK,a\l 



^K,aU = ^ (ti<x' - UK)y' 



1 - ■ {Xa - Xk) ] nK,a if CJ = d' 



\a'\nx ^1 ■ {Xfj — XK)nKa otherwise 



(28) 



Thus, 



/ Vvu{x) ■ A{x)Vvv{x)dx = J2 12 Yl ^if {u^ - UK){va' - vk) yu,veXv, (29) 



with, 



^k' = Yl • ^Kyy''"''' and Ak,." = I Aix)dx. (30) 

Then we get that the local matrices [A'^ )aa'&£K symmetric and positive, and the identification 
of the numerical fluxes using (21) leads to the expression: 

FkAu)= ^K'iuK-u,'). (31) 

Remark 2.8 (Link with the MFD method) The above technique yields an explicit construction 
of a particular MFD method. Indeed, if one chooses xk o,s the centre of mass of K , the matrix Ak 
defined by (30) is an adequate choice for the matrix We which is a parameter in the general formulation 
of the family of MFD methods as proposed in [10]. The advantages of the specific matrix Ak are that: 

• on particular meshes, taking a natural choice for xk (for instance the circumcenter for a tri- 
angular mesh in the case of a 2D isotropic problem), it degenerates to a diagonal matrix (see 
Lemma 2.1 below); 

• it is linked to an explicit formulation of a consistent gradient, which is used to define the discrete 
bilinear form (29) . 

Note however that the SUSHI scheme defined by (18) is not the MFD method of [9, 10]; the main 
reason is that according to the choice B, the SUSHI scheme may be either a completely cell-centred 
scheme, or a partly or fully hybrid scheme, while the MFD method is a pure hybrid scheme. Note 
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also that in SUSHI, one may take any point in cell K for x^, while the MFD schemes [9, 10] are 
constructed with the centre of mass (however, this choice might he generalised) . 

Note that the procedure which we describe in Section 2.2 to write a cell-centred scheme could he applied 
to any mimetic scheme (or low order mixed finite element scheme) to yield a centred scheme. However, 
further investigations are needed to determine under what conditions the present convergence analysis 
extends to mimetic schemes, and conversely, whether the mimetic analysis applies to the SUSHI scheme 
(ongoing work, [16]). 

The fluxes defined by (22)-(31) satisfy certain properties whicli are detailed in Lemma 4.4, and which 
allow us to prove the convergence of the scheme, as is shown in Theorem 4.1. Note that it seems difficult 
to deduce such properties from fluxes obtained by using natural expansions of regular functions. Note 
also that both Lemma 4.4 and Theorem 4.1 hold for general heterogeneous, anisotropic and possibly 
discontinuous fields A, for which the solution u of (6) is not in general more regular than u S Hq{Q). 
In the case where A and u are regular enough, the local flux consistency satisfied by (31) is used 
in order to obtain an error estimate, see Theorem 4.2. The coefficient \/d may be replaced by any 
positive real number without any change in the proof of convergence; in fact, for certain problems it 
can be interesting to use another coefficient, as described in [20] for the so called "SUSHI-P" scheme 
(P for parametric, meaning that the user may choose the stabilisation coefficient as well as the set of 
edges B). The choice Vd is however natural in the sense that with this value, if = 0, the scheme 
boils down in two dimensions to the well-known harmonic averaging five points scheme on rectangles 
and a four-point scheme on triangles; more generally, in any space dimension, even ii B ^ 9 and taking 
the most natural value for Uq- if cr E B, the resulting flux is a two-point flux on meshes that satisfy 
the "superadmissibility condition" (32), not necessarily with a harmonic averaging of A in the case 
a £ B; this is proven in the next lemma. Note that this superadmissibility condition is also satisfied 
by rectangular parallelepipeds in three dimensions but unfortunately not by tetrahedra. 

Lemma 2.1 (Superadmissible mesh and two-point flux) LetT> = {MjSjV) be a discretisation 
of Q in the sense of Definition 2.1, satisfying the following superadmissibility condition: 

UK a = ^^^^^ VA' G 7W, Va G £k. (32) 

dK,a 

Let us furthermore assume that A.{x) = A(cc)Id, where X is a piece-wise constant function from to 
M, which is equal to a constant Xk in each K € A4; then, the inner product defined by (21)-(25) reads: 

{U,v)p= ^ Xk -^^{UK - Ua){vK - Va). 

Moreover, choosing thanks to (32), x„ = {dK,aXL + dL^a^K) / {dK,a + c^l,(j) for a G <5int with M.^ = 
{K,L} in (15), the scheme (18) is the following two-point flux scheme: 

V FK,a = [ fix) dx, (33) 



K&M 



XxXlidx^a + dL,a) Cr \ -f ^ c 1-1 KA S T^ lr,A\ 

FK,a = — 3 ' , 3 —7 — [UK - ul) if a € tint riH, Ma = {K,L}, (34) 

^Kd-L,a + ^L(lK,a "X.cr + "L,tT 

j,^^^ ^ dK,.XK + dL,aXL W\ .^^^ ^^^^ nB, Ma = {K, L}, (35) 

aK,a + dL,a 0,K,a + "L.ct 

Fk,<t = Xk^^uk if (7 e £extr]£K, (36) 

aK,a 
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Proof. Let us compute (n, v)p under the assumptions of Lemma 2.1. From (21) and thanks to (27) 
we get: 



{u,v) 



V 



\] Xk / V©u(a;) • 'V-jyv{x)dx 



= ^ i \K\VkU ■ VkV + ^ RK,aU RK,aV j . 

Now from the definition (22) and thanks to the assumption (32), the discrete gradient given by (22) 
may be written as follows: 



From (23), we get 

>p 14 



E\(y\dK,cj ^fd , \ / \t Ir^lTJ 

— — ^- — [x„ - Xk) -: — [Xa - Xk) = \K\ld. 



Therefore, we get that 



^^^^^''^ RK,aU RK,aV = ^ ~ '"^^^'"'^ ~ ~ \K\\7kU ' ^RV, 



which in turn yields that 

{U,v)j,= ^ Xk ^ -j^iUa -UK){Va -Vk)- 

Hence the matrix Ak only contains the terms on the diagonal, and the fiux FK.aiu) is given by 

FK,a{u) = XK-^r^{uK - Ua). 

Then the scheme (18) can be written as a classical cell-centred finite volume scheme, with two-point 
fluxes Fk,l{u) = —Fk,l{u) for any a G <?int with Ai^ = {K, L]. Indeed, in the case a ^ B, the above 
expression of Fk,(t{u) allows us to get the following expression of Ua- from (9): 



^^^^^^^ 
+ 



This yields the harmonic averaging two-point flux 

7-1 / \ II fj dj^ (J / ^ 

= '^' ^ I ^ ^^^ " 

In the case a £ B, the two-point barycentric formula = {dK,aUL + dL,aUK)/idK,a + dL,a) together 
with (18) leads to the resulting two-point flux 

^^y^ _ dK,aXK + dL,aXL \<j\ _ 

dK,a + dL,a dK,c7 + C?L,o- 

□ 
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3 Numerical results 



We present some numerical results obtained with various choices of B in the scheme (18), (13) with 
the flux (31), which we synthesise here for the sake of clarity: 

Find u G X-p^^ (that is {uk)k&M': {ucr)cr<^n)j such that: 

FK,a{u){vK - Va) = ^ VK f{x)dx, for all V e Xt>,b, 

with FK,a{u) = ^ A^j^' {ua, -uk) yK e M,ya e £k- 

where the matrices are defined by (30)-(28). In the following, we shall use the choices i3 = 
(HFV), B = iSint or B the set of edges which are located on the diffusion tensor discontinuity interfaces; 
this latter choice is reported as SUSHI-NP (for non parametric) in [20], in contrast with SUSHI-P (for 
parametric) where the choice of the set B may be different, along with the value of the stabilisation 
coefficient in (25). 



3.1 Implementation 

Let us first describe an implementation aspects of the scheme. The unknowns, i.e. the values uk, 
iov K ^ M. and the values Uo-, cr £ £iat H are ordered as (Mi)i=i,...,Ar. The N x N matrix and the 
X 1 right-hand-side of the linear system resulting from (18) are computed thanks to a loop over 
the control volumes K & M and to an inner loop on each edge a G £k- Let us detail the matrix 
computation loop. 

1. All stored matrix coefficients are initially set to 0. 



2. The expression FK,a{u) is written in the form F^^aiu) = X]i=i N ^^K a'^i' where the nonzero 
coefficients {a^^ ^)i=i,...,N are only locally computed (they are not stored for all K and a). These 
coefficients are obtained after the elimination of all {ua)aeSKr]l3 iir (31): 



3. The line of the matrix corresponding to the unknown uk is incremented at the column j with 
the coefficient o!"^^. 

4. If (T G i3 with Va = YliLeM Pa'^L for any v G Xxt^B, the line of the matrix corresponding to each 
L ^ M such that 7^ is incremented at the column j with the coefficient —j3^a"^^. 

5. If o" G £\nt n the line of the matrix corresponding to the edge a is incremented at the column 
j with the coefficient —aj/^. 

This procedure is identical in the cases B = % (HFV), S 7^ and B = £\nt. However, in the case where 
S = (HFV), one may eliminate the unknowns uk with respect to the unknowns u^, as in the hybrid 
implementation of the mixed finite element method. 



3.2 Order of convergence 

We consider here the numerical resolution of Equation (1) supplemented by the homogeneous Dirichlet 
boundary condition (2); the right-hand side is chosen so as to obtain an exact solution to the problem 
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and easily compute the error between the exact and approximate solutions. We consider Problem 
(l)-(2) with a constant matrix A: 

= ( il ) • P«) 

and choose / : $7 — > M such that the exact solution to Problem (l)-(2) is u defined by u{x,y) = 
16x(l — x)y{l — y) for any (rc, y) G (7. Note that in this case, the composite scheme is in fact the 
cell-centred scheme, there are no edge-unknowns. 

Let us first consider conforming meshes, such as the triangular meshes which are depicted on Figure 
1, and uniform square meshes. 



Figure 1: Regular conforming coarse and fine triangular grids 

For both B = 9 (pure hybrid scheme: HFV) and B = £i^t (cell-centred scheme), the order of conver- 
gence is close to 2 for the unknown u and 1 for its gradient. Of course, the hybrid scheme is almost 
three times more costly in terms of number of unknowns than the cell-centred scheme for a given 
precision. However, the number of nonzero terms in the matrix is, again for a given precision on 
the approximate solution, larger for the cell-centred scheme than for the hybrid scheme. Hence the 
number of unknowns is probably not a sufficient criterion for assessing the cost of the scheme. 
Results were also obtained in the case of uniform square or rectangular meshes. They show a better 
rate of convergence of the gradient (order 2 in the case ofTi = Si^t and 1.5 in the case B = £int), even 
though the rate of convergence of the approximate solution remains unchanged and close to 2. 
We then use a rectangular nonconforming mesh, obtained by cutting vertically the domain into two 
parts and using a rectangular grid of 3n x 2n (resp. 5n x 2n) on the first (resp. second side), where 
n is the number of the mesh, n = 1, . . . , 7. Again, the order of convergence which we obtain is 2 for u 
and around 1.8 for the gradient. We give in Table 1 below the errors obtained in the discrete norm 
for u and Vti for a nonconforming mesh and (in terms of number of unknowns) and for the rectangular 
4x6 and 4 x 10 conforming rectangular meshes, for both the hybrid and cell-centred schemes. We 
show in Figure 2 the solutions for the corresponding grids (which look much the same for the two 
schemes) . 

Further detailed results on several problems and conforming, nonconforming and distorted meshes 
may be found in [20]. 

3.3 The case of a highly heterogeneous tilted barrier 

We now turn to the heterogeneous case. The domain =]0, l[x]0, 1[ is composed of 3 sub-domains, 
which are depicted in Figure 3: Jli = {{x,y) £ fl; ipi{x,y) < 0}, with ipi{x,y) = y — 5{x — .5) — .475, 
^^2 = {{x,y) G n;ipi{x,y) > 0,ip2{x,y) < 0}, with ip2{x,y) = ipi{x,y) -0.05, ^3 = {{x,y) G 
il; (p2{x, y) > 0}, and 5 = 0.2 is the slope of the drain (see Figure 3). Dirichlet boundary conditions are 
imposed by setting the boundary values to those of the analytical solution given by u{x, y) = —ipi (x, y) 
on U 1^3 and u{x,y) = — (^1(2;, y)/10~^ on $72- 
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NU 


NM 


eiu) 


e(V'u) 


n 


S = 13 = £int 


S = I3 = £int 


B = % B = Sint 


6 = 13 = £int 


CI 

NC 
C2 


130 48 
182 64 
222 80 


874 488 
1334 724 
1542 864 


1.28E-01 1.20E-01 
1.03E-01 9.43E-02 
7.61E-02 7.09E-02 


1.64E-02 3.57E-02 
1.66E-02 3.69E-02 
9.18E-03 2.44E-02 



Table 1: Error for the nonconforming rectangular mesh, pure hybrid scheme {B = 0) and centred {B ~ £int) 
schemes. For both schemes NU is the number of unknowns in the resulting linear system, NM is the number of 
nonzero terms in the matrix, e(u) is the discrete norm of the error of the solution and e(VM) is the discrete 
norm of the error in the gradient. CI and C2 are the two conforming meshes represented on the left and 
the right in Figure 2, and NC is the nonconforming one represented in the middle. 




Figure 2: The approximate solution for conforming and nonconforming meshes. 
8x6 mesh, centre: nonconforming 4 x 6, 4 x 10 mesh, right: conforming 10 x 10. 



Left: conforming 



The permeability tensor A is heterogeneous and isotropic, given by A{x) = A(a;)Id, with X{x) = 1 for 
a.e. X G ill U ils and X{x) = for a.e. x £ 0.2- Note that the isolines of the exact solution are 
parallel to the boundaries of the sub-domain, and that the tangential component of the gradient is 
0. We use the meshes depicted in Figure 3. Mesh 3 (containing 10 x 25 control volumes) is obtained 
from Mesh 1 by the addition of two layers of very thin control volumes around each of the two lines 
of discontinuity of A: because of the very low thickness of these layers, equal to 1/10000, the picture 
representing Mesh 3 is not different from that of Mesh 1. 




Figure 3: Domain and meshes used for the tilted barrier test: mesh 1 (10 x 21 centre), mesh 2 (10 x 100 
right) 

We get the following results for the approximations of the four fluxes at the boundary. 
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nb. unknowns 


matrix size 


X = 


X = 1 


y = 


y = l 


analytical 






-0.2 


0.2 


1. 


-1. 




mesh 1 


210 


2424 


-1.17 


1.17 


3.51 


-3.51 


B = Sint 


mesh 2 


1000 


11904 


-0.237 


0.237 


1.104 


-1.104 




mesh 3 


250 


2904 


-0.208 


0.208 


1.02 


-1.02 


SUSHI-NP 


mesh 1 
mesh 2 


239 
1020 


2583 
12036 


-0.2 
-0.2 


0.2 
0.2 


1. 
1. 


-1. 
-1. 


HFV 


mesh 1 


599 


4311 


-0.2 


0.2 


1. 


-1. 


mesh 2 


2890 


21138 


-0.2 


0.2 


1. 


-1. 



Note that the values of the numerical solution given by the pure hybrid (HFV) and composite (SUSHI- 
NP) schemes are equal to those of the analytical solution (this holds under the only condition that 
the interfaces located on the lines ifi{x, y) = 0, i = 1, 2, are not included in B, and that, for all a £ B, 
all K e M with /3f / are included in the same sub-domain J7j). Note that Mesh 3, which leads 
to acceptable results for the computation of the fluxes, is not well suited for such a coupled problem, 
because of too small control volume measures. Hence SUSHI on Mesh 1 appears to be the most 
suitable method for this problem. 

A satisfying natural choice (SUSHI-NP in the above results) is thus to match 7i with the discontinuities 
of A. It is sometimes interesting to choose another set B. This is for instance the case for the numerical 
locking problem for which the choice ;B = is best even though the diffusion tensor is homogeneous 
[20]. 

It is also sometimes interesting to replace the stabilisation coefficient Vd in (25) by some other coef- 
ficient a > 0. This is the case for instance for very distorted meshes or singular problems, in order 
to maintain the positivity of the unknown. The coefficient a is taken to be greater than ^/d. The 
approximate solution remains positive, but the norm of the error is generally larger. We refer to 
[20] for such experiments. 



4 Convergence of the scheme 

Let us first introduce some notations related to the mesh. Let P = {M,£,V) be a discretisation of Q 
in the sense of Definition 2.1. The size of the discretisation D is defined by: 

h-p = sup{/ix, K E M}, 

and the regularity of the mesh by: 

9d = max ( max , max | • (39) 

For a given set B C <fint and for a given family {[3^) kgm satisfying property (15), we introduce a 
measure of the resulting regularity by 

Ot> B = max 0x>, max — . (40) 

\ K£M,aGeKnB I 

Remark 4.1 Note that, for any mesh, it is easy to choose the family ij3^) ksm so that 6x),b remains 

small. It suffices to express x„ as the bary centre of d + 1 points xl (which is always possible), for L 
sufficiently close to K , so that — Xfj is close to Hk when (3^ ^ 0. Note also that in fact, it would 
be sufficient to have h\ with rj > 1 instead of h\ in (40) thus allowing the use of farther points. 
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Remark that, thanks to the assumption that K is a; /^-star-shaped, the following property holds: 

^ \(T\dK,a = d \K\ ^KeM. (41) 

The space Xx> defined in (10) is equipped with the following semi-norm: 

yvGXj,, \v\\= E j^K-^^)'' (42) 

which is a norm on the spaces X-ofi and Xt>^b respectively defined by (11) and (17). 
Let Hm{^) C L2(0) be the set of piece-wise constant functions on the control volumes of the mesh 
M. We then denote, for all v G Hf^{Q) and for all a G with ^A(J = {K, L}, D^jV = \vk — vl\ and 
d(T = dK,a + c^L,cr, and for all a G Sext with A4a = {K}, we denote D(^v = \vk\ and d^ = dK,a- We 
then define the following norm: 

^vGHMi^), \\vh,2M= E E Hd/^,.(^)' = El^l^^- (43) 
(Note that this norm is also defined by (74) in Lemma 5.2, setting p = 2). 

For all V £ X^i, we denote by Hm^ ^ the piece- wise function from i7 to M defined by 

YlM^ix) = vk for a.e. x G K, for all K G M. Using the Cauchy-Schwarz inequality, we have for all 
a G fint with Ma = {K, L}, 

1 < \ V-u G Ax), 

da dK,a dL,a 

which leads to the relation 

\\^Mv\\l2,M < Hi G Xt>,o. (44) 

For all ip G C(r2,M), we denote by Pxnp the element of Xjy defined by {{ip{xx))KeM^ {^{XcT))aG£)i 
by Pv,bV the element v G X-^^b such that vk = ^{xk) for all K G A^, Uo- = for all a G fextj 
Va = YIki^M f^a^^i^K) fo^ all fj G ;B and Va = ^{Xa) for all o" G 

We denote by PmV ^ Hm{Q) the function such that ^^(/^(a;) = if^xx) for a.e. a; G .fC, for all K £ 
(we then have PmV' = ^M^v^ = n_A4PD,BV9). 

The following lemma provides an equivalence property between the L^-norm of the discrete gradient, 
defined by (22)- (26) and the norm | • |x- 

Lemma 4.1 Let D be a discretisation of Q in the sense of Definition 2.1, and let 9 > 9x> be given 
(where Ojy is defined by (39) J. Then there exists Ci > and C2 > only depending on 9 and d such 
that: 

Ci < ||VDn||i2(Q) < C2 |n|x G Xx,, (45) 

where Vxi is defined by (22) -(26). 



Proof. By definition, 



Therefore, using property (27), 



l|Vl.u|li.(^). = E ( \K\\^KU\' + E ""^^{RK^auf ) . (46) 



K&M \ a££K 
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Let us now notice that the foHowing inequahty holds: 

A 



{a - bf > Y^a^ - Xb"^ Va, 6 G M, VA > -1. (47) 



We apply this inequality to {RK,au)'^ for some A > and obtain 
This leads to 

.fee ^ ^ + ^.fe ^ ^^"^ ^ 

Choosing A = we get that 

which shows the left inequality of (45). 

Let us now prove the right inequality. On one hand, using the definition (22) of V/^n and (41), the 
Cauchy-Schwarz inequality leads to 

K O-K.a \J^\ d'K,cr 

On the other hand, by definition (25), and thanks to the definition of the regularity of the mesh (39), 
we have 

{RK,^uf < 2d {i^^^^f + l^^^l'l^^^l') ^ 2d ((^^^^)' + ^'iVi^np) . (50) 

From (46), (49) and (50), we conclude that the right inequality of (45) holds. □ 
We may now state a weak compactness result for the discrete gradient. 

Lemma 4.2 (Weak discrete compactness) Let he a family of discretisations in the sense 
of Definition 2.1 such that there exists > with 9 > 9x> for all D G JT. Let {ut>)v(^j^ be a family of 
functions, such that: 

• ux> £ Xx>fi for all V £ 

• there exists C > with \ud\x < C for all V £ 

• there exists u £ L'^{fl) with lim ||n_A4UD — ullLacf^-) = 0. 

Then, u £ H^{Vl) and Vx>ux> weakly converge in L'^{Q)'^ to Vti as hx> — > 0, where the operator Vxi is 
defined by (22) -(26). 

Proof. Let us prolong Hj^ud and Vdud by outside of fi. Thanks to Lemma 4.1, up to a 
subsequence, there exists some function G £ L'^iW^)'^ such that \7x>ux> weakly converges in L^(M'^)'^ 
to G as hv ^ 0. Let us show that G = Vu. Let i/> £ C^{W^)'^ be given. Let us consider the term 
defined by 

= / Vx>uv{x) ■'ip{x)dx. 
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We get that Tf = Ti^ + T^, with 

= ^ \cr\{u^ - UK)nK,a ■ '4>K, with iI)k = / ip{x)d. 

and 

We compare T2 with T^^ defined by 



with 



We get that 



K&M o&Sk 

■0<T = rr / •j/'(^c)d7(a;). 



which leads to hm {T^ - rf ) = 0. 
Since 

K&M a££K 

we get that hm = — u{x)divrp{x)dx. Let us now turn to the study of T^. Noting again that 
(27) holds, we have: 

Since ■?/> is a regular function, there exists only depending on ■0 such that | J^^^ iip{x) — il}K)dx\ < 
Cxjjhxr — From (50) and the Cauchy-Schwarz inequality, we thus get: 

lim T? = 0. 

This proves that the function G E L^(M'^)'^ is a.e. equal to Vu in M'^. Since u = Q outside of 0, 
we get that u G Hq{Q), and the uniqueness of the limit implies that the whole family Vvut> weakly 
converges in L'^iW^)'^ to Vn as /i© 0. 
□ 

Note that the proof that u G HQiyi) also results from (44), which allows us to apply Lemma 5.7 of 
the Appendix in the particular case p = 2. Let us also remark that several discrete gradients could 
be chosen, which satisfy the weak compactness property (see for instance the proof of Lemma 5.7). 
However, we emphasise that the choice of the specific gradient (22) also stems from coercivity and 
consistency issues. Let us now state the discrete gradient consistency property. 

Lemma 4.3 (Discrete gradient consistency) Let D be a discretisation of Q in the sense of Defi- 
nition 2.1, and let 9 > 9x> be given. Then, for any function ip G (7^(0), there exists C3 only depending 
on d, 9 and ip such that: 

- ^v\\{L^{n)Y < C'3 h-D, (51) 

where Vd is defined by (22) -(26). 
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Proof. From definitions (26) and (24) we get 
From (22), we have, for any K £ M, 

= j-^ ^ kl {V(p{xk) ■ {Xa - Xk) + hj^pK,a)nK,a, 

where 1/5^,0-1 < with only depending on ip. Thanks to (23) and to the regularity of the mesh, 
we get 

From this last inequality, using Definition 25, we get 

\RK,aPv'^\ 



which concludes the proof. □ 
We now give the abstract properties of the discrete fluxes, which are necessary to prove the convergence 
of the general scheme (18), (13), and then prove that the fluxes that we constructed in Section 2.4 
indeed satisfy these properties. 

Definition 4.1 (Continuous, coercive, consistent and symmetric families of fluxes) 

Let T he a family of discretisations in the sense of definition 2.1. ForT) = {M.,£,V) £ J-', K £ M and 
a € £, we denote by F^ ^ a linear mapping from X-jy to M, and we denote by ^ = {{F^ ^) KeM)v&r ■ 

We consider the bilinear form defined by 

{u,v)f = PkA^)^^k - ^("'^) e ^l- (52) 

The family of numerical fluxes $ is said to be continuous if there exists M > such that 

{u,v)f < M\u\x\v\x y{u,v) £ Xl, yv = {M,£,V) £ (53) 

The family of numerical fluxes ^ is said to be coercive if there exists a > such that 

a\u\j. < {u,u)f 'iuG Xv^V = {M,£,V) £ J". (54) 

The family of numerical fluxes $ is said to be consistent (with Problem (l)-(2)j if for any family 
{uv)v(iT satisfying: 

• uj) £ Xx>fl for all V £ J^, 

• there exists C > with ^ C for all T> £ 

• there exists u £ L'^{Q) with lim ||n_A4Ux) — ^illL^rrj) = (recall that, from Lemma 5.7, we get 
that u £ H^i^)), 



^ \ f{x^) - (p{xk) - VkPv^ ■ {xa - Xk) 



< ^{hj<PK,a + hl dC^e) 

< Vde{hKC^ + hKdC^9), 
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then 



lim {uv,Pvv)f = / A{x)Vip{x) ■ Vu{x)dx Vv? G C~(0). (55) 

Finally the family of numerical fluxes <I> is said to be symmetric if 

{u,v)f = {v,u)f y{u,v) G Xl, yv= {M,£,V) G J". 

We now show that the family of fluxes defined by (28)-(31) satisfies the definition of a consistent, 
coercive and symmetric family of fiuxes. Recall that the SUSHI scheme (37) is studied numerically in 
Section 3 with this choice for the family of fluxes. 

Lemma 4.4 (Flux properties) Let !F be a family of discretisations in the sense of Definition 2.1. 
We assume that there exists ^ > with 

e^<e yv = {M,£,v) e (56) 

where 9x> is defined by (39). Let $ = {{F'^ ^) kgm)v£T be the family of fiuxes defined by (28)-(31). 

Then, the family ^ is a continuous, coercive, consistent and symmetric family of numerical fluxes in 
the sense of Definition 4-F 

Proof. Since the family of fiuxes is defined by (28)-(31), it satisfies (21), and therefore we have: 

{u,v)f= / Vvuix) ■ A{x)'VT>v{x)dx \/u,v £ Xx). 
Jn 

Hence the property {u,v)f = {v,u)f holds. The continuity and coercivity of the family $ result from 
Lemma 4.1 and the properties of A, which give: {u,v)f < -^||Vx>ti||2,2(Q) || Vx)w||j;^2(f^) and {u,u)f > 
A|| Vx'ii||^2(f^) for any u,v £ Xx>. The consistency results from the weak and strong convergence 
properties in Lemmas 4.2 and 4.3, which give \/x>uv — > weakly in L'^{Q) and \/j)Px>(p V(/9 in 
L^(r2) as the mesh size tends to 0. □ 

Theorem 4.1 (Convergence) Let T be a family of discretisations in the sense of Definition 2.1, 
for any D £ T , let B C Suit o,rid {(3^) k^m satisfying (15). Assume that there exists > such that 

Gv,B ^ ^) for all D £ T , where 9t>,]3 is defined by (40). Let # = {{F^ ^)KeM)'D<^T be a continuous, 
coercive and symmetric and consistent family of numerical fluxes in the sense of Definition 4-1- Let 
{ux>)t>gj^ be the family of functions satisfying (18) for all V £ . Then H^ux) converges in L'^{Q) to 
the unique solution u of (6) as hv 0. Moreover 'Vx>ut) converges to Vu in L'^{VtY as hv 0. 

Proof. Letting t; = in (18) and applying the Cauchy-Schwarz inequality yields 

{uv,uv)f= / f{x)UMUvix)dx < \\f\\L2tQ)\\IlMUv\\L^n)- 
Jn 

We apply the Sobolev inequality (77) with p = 2, which gives in this case 

l|nAinD||j;^2(Q) < C4||n_A4u||i,2,A^. 

Using (44) and the consistency of the family $ of fluxes, we then have 

a\UMUv\x < C'4 ll/llL2{n)I^^X)U- 

This leads to the inequality 

\\uv\\i,2,M < Wv\x < -^Wfh^in)- (57) 
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Thanks to Lemma 5.7, we get the existence of n G ffQ(r2), and of a subfamily extracted from T ^ such 
that ||n_A^UD — ^11^2(^2) tends to as /i© — > 0. For a given ip G C^(0), let us take v = Pv,B^ in (18) 
(recall that Px> ^ ^v,b)- We get 



{uv,Pv,b^)f = / f{x)PM^{x)dx. 
Jn 



Let us remark that, thanks to the continuity of the family ^ of fluxes, we have 

{uv,Pv,B^ - Pvip)F < Af— — ||/||i2(Q) \Pv,B^ - Pv^\x- 

Thanks to (15) and (40), we get the existence of only depending on ip (through its second order 
partial derivatives) such that, for all -fT G and all o" G ;B H £k, 

I ^ Pj^ifixL) - ip{x„)\ < \Po\\xL - x„?C^ < ev,BC^h\. (58) 
We can then deduce 

lim \Pv,BV-Pv^\x = 0. (59) 

Thanks to the .7-"-extracted subfamily properties, we may apply the consistency hypothesis on the 
family <I> of fluxes, which gives 



h-p 

Gathering the two results above leads to 



lim {ut>, Px>p)f = / A{x)'Vip{x) ■ 'Vu{x)dx. 
■v->o Jn 



which concludes the proof of the following equality 



Yim. {uv , Pv ,bv) F = I h.{x)VLp{x) ■Vu{x)dx, 



K{x)V (p{x) ■ Vu{x)dx = I f{x)ip{x)dx. 

Jn 



Therefore, u is the unique solution of (6), and we get that the whole family {uT>)veJ^ converges to u 
as hi) — > 0. 

Let us now prove the second part of the theorem. 

Let ip G C^(il) be given (this function is devoted to approximate u in i?^(0)). Thanks to the 
Cauchy-Schwarz inequality, we have 



/ \Vvuv{x) - Vu{x)\^dx < 3 (T^ + Ti" + Tj), 
Jn 



withT^ = J^\Vvuv{x)-VvPvv{x)\^dx,T^ = J^\VvPvV^{x) -Vpix)\^dx, and T7 = J^\Vpix) 
S/u{x)\'^dx. Thanks to Lemma 4.3, we have liuih^^oT^ = 0. 

Thanks to Lemma 4.1 and to the coercivity of the family of fluxes, there exists C5 such that 

||Vx,v||^2(Q)d < C2^\v\x < C^{v,v)f Vt; G Xv, 
with C5 = Taking v = ux> — Pv'P, we have 

< C5 {{uv,ux>)f - 2(ux), Pvv)f + {Pv^, Pv^)f)- 
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By Theorem 4.1 and thanks to and consistency of the family of fluxes, we get 

lim {ux>,Pv'p)f = / 'Vu{x) ■ A{x)'Vip{x)dx. 



The sequence \Pv^\x is bounded; using the regularity of ip, the regularity hypotheses of the family of 
discretisations, together with the consistency of the family of fluxes implies that 



lim {Pd if, Pv'p)f = / Vifix) ■ A{x)'V(p{x)dx. 



Remarking that passing to the limit hj) ^ in (18) with v = ut> provides that {uv,ux>)f converges 
to Vit • AVudx, we get that 

lim {uj} — Pvf, ux> — Pvf) F = V('u — (p) ■ AV(n — ip)dx < A / | Vn — V^^^da;, 
which yields 

limsupTg^ < Cs A / \Vu- V^pl'^dx. 
h-D-^o Jn 

From the above results, we obtain that there exists Cg, independent of T>, such that 



/ \Vvuv{x) - Vu{x)fdx < Ce / iV^pix) - Vu{x)\'^dx + 
Jn Jn 



8 ; 



with (noting that f is fixed) lim/j^^o = 0. Let e > 0; we may choose ip such that \ Vip{x) — 
Vu{x)\'^dx < e, and we may then choose /i© small enough so that < e. This completes the proof 
that 

lim / \Vx>uv(x)-Vu(x)\^dx = (60) 

hT,^0Jn 

in the case of a general continuous, coercive, consistent and symmetric family of fluxes. □ 

Let us write an error estimate in the particular case A = Id, assuming a regular exact solution to (6). 

Theorem 4.2 (Error estimate, isotropic case) We consider the particular case A = Id, and we 
assume that the solution u G Hq{P,) of (6) is in (7^(0). Let V> = {M,£,V) be a discretisation in the 
sense of Definition 2.1, let B C £int be given, let B = {l3^)rj^B,K&M C K such that (15) holds, and let 
S ^ (^v,B be given (see (40)j. Let {FK,a)KeM,ae£ be a family of linear mappings from Xx> to M, such 
that there exists a > with 

(Av\'x < {v,v) F E Xx), (61) 

defining {■,-)f by (52). We denote by 

^(^)= E E ^ f ^A>(Pp,en) + / Vu{x)-nK^Alix)\ . (62) 

Then the solution ux> of (18) satisfies that there exists C^, only depends on a and on 6, such that 

W^MUv - PMu\\L-i(n) < CtE{u), (63) 
and satisfies that there exists Cg, only depending on a, 6 and u such that 

\\Vvuv - Vn||i2(o)<i < {E{u) + hv) ■ (64) 

Moreover, in the particular case where {FK,a-)K&M,<T&£ is defined by (28)-(31), there exists Cg, only 
depending on a, 9 and u, such that 

E{u)<C<shv. (65) 
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Remark 4.2 (Extensions of the error estimate) Note also that the extension of Theorem 4-2 to 
the case u £ H'^{^) is possible for d = 2 or d = 3. However it would demand a rather longer and 
more technical proof and is not expected to provide more information on the link between accuracy 
and the regularity of the mesh than the result presented here. In the case of the pure hybrid scheme 
(HFV, B = $), an error estimate could however be obtained by assuming u piece-wise to be H^. 
Such error estimates were also obtained for pure hybrid schemes of the mimetic type by using the 
tools of the mixed finite element theory (see e.g. [10]). If B ^ ^, one must furthermore assume 
that the barycentric formulae (14) -(15) or (19) -(20) are written with unknowns located in the same 
regularity zone, as explained in Remark 2.6. Nevertheless such error estimates are not possible for 
general L°° diffusion operators, since in such a case the maximal regularity of the continuous solution 
is Hq{Q). Then, by interpolation, one may get some error estimates if the continuous solution is in 
IIq{Q) nH^{Q) as in the classical finite element framework. 

Proof. Let v G Xxi, since —Au = f, we get: 

-5:„./A„Wd.= //(.)n^.Wd.. (66) 

Thanks to the following equality (recall that u G C^(17) and therefore Vu • nK,a is defined on each 
edge a) 

- ^ j Au{x)diX = - ^ "^{vK-Va) / Vu{x) ■nK,Al{x), 

we get that 

{Pv,bu,v)f= / f{x)UMv{x)dx + ^ ^ (f^^^{Pv,bu) + / \/u{x) ■ nK,ad'y{x)\ {vk - v^). 
Taking v = Pt>,bu — n-p G in this latter equality and using (66) we get 

{v,v)f= ^ J2 (^^,^7(^1',^^)+ / "^Uix) ■nK,ad'y{x)\ (VK -Va), 

which leads, using (61) and the Cauchy-Schwarz inequality, to 

a\v\x < E{u). (67) 

Using (44) and the Sobolev inequality (77) with p = 2 provides the conclusion of (63). Let us now 
prove (64). We have 

\\Vvuv - Vn||i2(Q)d < \\Vvuv - Vx)-Px),Bn||i2(t^)d + \\^vPv,bu - Vti||^2(t^)d. 

The bound of the first term in the above right-hand side is bounded thanks to Lemma 4.1 and (67). 
The inequality \\^vPv,BU — Vti||j;^2(Q)<i < Ciq/i© is obtained thanks to Lemma 4.3 and using a similar 
inequality to (58), replacing ip by u. 

Let us now turn to the proof of (65) in the particular case where the family of fluxes is defined by 
(28)-(31). Indeed, we get in this case that, for all v G Xj), 
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with 



y 



li^l 
|cr| 



li^l 



L \K\ dK,a'\K\ 



otherwise . 



Using (23), we get that 



2^ j^y'''' = \a\nK,a. 



Since there exists Cn £ M+ such that |-R_ft-,a-'-pD,B^^| < Cu hx, there exists some C12 £ IR+ with 



FK,a(.Pv,Bu) + / Vn(a;) • nj^^^d7(a;) 

J a 

This leads to the conclusion of (65). □ 



< C12 \(T\hK- 



5 Discrete functional analysis 

This section is devoted to some results of functional analysis that are useful for the proof of convergence 
of numerical schemes when the approximate solution is piece-wise constant on the mesh. Although 
some of the results presented here were already introduced in previous works of the authors, they 
were mostly presented (even when not needed, see [17, Remark 9.13 p. 793]) in the framework of 
"admissible" meshes, that is meshes with an orthogonality condition. 

We recall that in the proof of the main convergence Theorem 4.1, we first obtain from the scheme 
some estimates on the approximate solutions in the discrete norm. We now show how, from a 
general discrete P^^'P estimate (this generalisation to p 7^ 2 is useful in the case of nonlinear problems) 
we obtain a discrete L'^ estimate for some q > p (Lemma 5.3). We then obtain a certain compactness 
result in (Lemma 5.5 and therefore in (Lemma 5.6), which in turn allows to show that the limit 
of the approximate solution is in VFQ'''(r2) (Lemma 5.7). 



5.1 Discrete Sobolev embeddings 
5.1.1 Discrete embedding of W^'^ in 

The discrete Sobolev embedding of W^'^ in requires less assumptions on the mesh than those 
given in Definition 2.1. We therefore introduce a larger class of meshes in the following definition. 

Definition 5.1 (Polyhedral partition of il) Let d > 1 and let be an open bounded set in W^, 
whose boundary is a finite union of part of hyperplanes. A polyhedral partition Ai of Q is a finite 
partition of ft such that each element K of this partition is measurable and has a boundary dK that 
is composed of a finite union of parts of hyperplanes (the facets of K) denoted by a: dK = Uo-ef^^ cr. 
Let £ be the set of the facets of all the elements of M: £ = Ukgm^k ■ If (y ^ £ is a facet of this 
partition, one denotes by \cr\ the {d — 1)-Lebesgue measure of a. Let iif^(O) be the set of functions 
from to M, constant on each element of M. Let u S II m{^)- If ^ H £l (that is a is a facet 
such that a C K n L), one sets D^u = \uk — ul\- If a ^ £ is on the boundary ofQ, and K £ M. (that 
is a = r) K), one sets D„u = \uk\- For u G II_m{0,), one sets 

\\u\\i,i,M = '^W\F>aU. (68) 
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Lemma 5.1 Let d > 1 and let Q be an open bounded set ofM.'^, whose boundary is a finite union of 
parts of hyperplanes. Let M be a polyhedral partition of Q in the sense of Definition 5.1. Then, with 
the notations of Definition 5.1, 

||w|lLi*(n) < ^^ll^lli'i'-^ ^ Hm{^), with 1* = -j^- (69) 



Proof. Different proofs of this lemma are possible. A first proof consists in adapting to this discrete 
setting the classical proof of the Sobolev embedding due to L. Nirenberg (actually, it gives 1/2 instead 
of l/{2yjd) in (69)): it is based on an induction on d. This proof is essentially given in [17, Lemma 9.5 
page 790], with slightly less general hypotheses; in fact the so called orthogonality assumption is not 
used in the proof of Lemma 9.5 of [17]. An easy adaptation of this proof leads to the present lemma 
(with 1/2 instead of l/{2y/d) in (69)). 

The present proof makes direct use of L. Nirenberg's result, namely: 



< 7T^||^|Ih^M(m<^) Vtx G W^'\W'), (70) 



1 

2d' 

where ll""!! vk1'1(R'*) — SiLi ll-^i^llLi(R'*) and DiU is the weak derivative (or derivative in the sense of 
distributions) of u in the direction Xi (with x = (xi, . . . , Xd) G W^). 

For u G L^iW^), one sets = with, for i = 1, . . . , d, = sup{/ n^da;, 

if G C^iW^), llvllLoo(Md) ^ !}• The function u belongs to the space BV\iu(^ L^{W^) and ||n||By < oo. 
We first remark that (70) is true with instead of ||ti||i4/i,i(Md), and \iu G instead of VF"'^'^(M'^). 

Indeed, to prove this result (which is classical), let p G C^(M'^,M_|_) with J pdx = 1. For n G N*, 
define p„ = n'^p{n-). Let u G BV and Un = u-k pn so that, with (70): 

1 

\\Un\\L^*(M.'') ^ -^'^\\DiUn\\L^(K'i)- (71) 
1=1 



Since Un is regular, ||i^iWn||Li(Rd) = ||Z)iM,i||M) and, for ip G C^(M ), using Fubini's theorem: 

f dip f d 

/ UnT—dx= / U——{ip-kpn)dx<\\Diu\\M\\^\\L°°(Rd). 

J^d dxi J^d dxi ^ ' 

This leads to ||L'jitn||ii(iRd) < [[DjuUM- Since Un u a.e., as n — > cx), at least for a sub-sequence, 
Fatou's lemma gives, from (71): 

„ „ 1 „ „ 

\\u\\Li*(Rd) < -^\M\bv G BV. (72) 

Let u G Hm{^)- One sets n = outside Vl so that u G L^(M'^). One has ||tt||sy = supj/jg^j u d\v(p dec, 
ip G C^(R"',M'^), ||v?||l°°{r<*) < 1}. with ||(/5||Loo(iRd) = supj^i^...^^ \Wi\\L^(wd) and p = {pi,.. .,pd)- But, 
for p G C^(M'^,M'^) such that ||¥'||_Lt=o(iRd) < 1, an integration by parts on each element of Ai gives 
(where is a normal vector to a and 7 is the {d — 1)— Lebesgue measure on a): 



I u diYp dx = 2_] D(tU / \p • T^aldjix) < Vd\\ 



^lll,l,A^- 



Then, one has ||u||_By < V^||^i||i,i,A^ and (72) leads to (69). 
□ 
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5.1.2 Discrete embedding of W^''p in L^* , I < p < d 

We now prove a discrete Sobolev embedding for 1 < p < d and for meshes in the sense of Definition 2.1. 

Lemma 5.2 Let d > 1, 1 < p < d and let Q. he a polyhedral open hounded connected suhset ofM.'^. Let 
D he a discretization on Q in the sense of Definition 2.1. Let r] > be such that rj < dK,a/dL,a < 1/?? 
for all a £ £, where = {K,L}. Then, there exists C13, only depending on d, p and rj such that 



< C'lS Ihlli.p.A^ Vn G Ht>{0.), (73) 



where p* = and 



\<.M = E El-I^A>f^)', (74) 



with da = dK,a + dL,a, if Ma = {K, L} , and d^ = dK,a, if Ma = {K}. 

Proof. We follow here L. Nirenberg's proof of the Sobolev embedding. Let a be such that al* = p* 
(that is a = p{d — l)/{d — p) > 1). Let u S Hx>{^). Inequality (69) applied with instead of u 
leads to: 

— 

( / \ufAx] ' <Y^\cj\Da\u\'' . 

For a G £i,,t, Ma = {K,L}, one has Dalul"' < a{\uK\°'~'^ + \uL\'^-'^)DaU. For a G fext, Ma = {K}, 
one has Da\u\°' < a\uK\°'~^L>aU. This yields: 

d-l 

I \ufdx\ ' < E E W\a\uKr~^DaU, (75) 



For ah a £ £, one has 1 < if a G Sint, Ma = {K,L}, or if cr G £c^t, Ma = {K}. Then, 

Holder's inequality applied to (75) yields, with q = p/{p — I). 

Since (a — l)q = p*, one has: 

E W\dK,a\nK\'^''-^^'' = d I \ufAx. 

Then, noticing that {d- l)/d- 1/q = 1/p*, we deduce (73) follows from (76) with C13 = a^^d^/i 
only depending on d, p and r]. □ 

5.1.3 Discrete embedding of W^'^ in L^', for some q > p 

Let 1 < p < 00, we now deduce from Lemma 5.3 the following lemma, which gives the discrete 
embedding of W^'^ in L'^, for some q > p. 

Lemma 5.3 Let d > 1, 1 < p < 00 and let Q he a polyhedral open bounded connected subset ofM.'^. 
Let T> he a mesh of in the sense of Definition 2.1. Let rj > be such that rj < dx a/dL,a ^ ^/v for 
all a G £, where Ma = {L^,L}. Then, there exists q > p only depending on p and there exists C14, 
only depending on d, 17, p and rj such that 

<Cia\\u\\i^p^m \/u£HD{n), (77) 

where \\u\\^ is defined in (74). 
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Proof. If p = 1, one takes q = 1* and the result follows from Lemma 5.1 (in this case C14 does not 
depend on t]). If 1 < p < d, one takes q = p* applies Lemma 5.2. 

I{ p > d, one chooses any q £]p,oo[ and pi < d such that p\ = q (this is possible since p\ tends to 
00 as pi tends to d). Lemma 5.2 gives, for some C13 only depending on p, d and r/, that ||ii||i9(n) ^ 
C13 ||u||i,pi,Ai- But, using Holder's inequality, there exists C15, only depending on d, p, O, such that 
||^^||i,pi,A^ ^ Ci5 \\u\\i^p^M- Inequality (5.3) follows with C14 = C13 C15 . □ 

5.2 Compactness results for bounded families in the discrete W^'^ norm 
5.2.1 Compactness in 

We prove in this section that bounded families in the discrete W^'^ norms are relatively compact in 
U'. We begin here also with the case p = 1, giving in this case a crucial inequality which holds for 
general polyhedral partitions of $7. 

Lemma 5.4 Let d > I and let Q, be an open bounded set in M.'^, whose boundary is a finite union of 
parts of hyperplanes. Let M be a polyhedral partition of VL in the sense of Definition 5.1. Then, with 
the notations of Definition 5.1, 

\\u{- + y)- u\\Li(ud) < \y\Vd\\u\\i^i^M Vn G Hm{^), Vy G R'^, (78) 



where u is defined on the whole space W^, taking u = outside Q., and \h\ is the Euclidean norm of 
h G M.'^. 

Proof. One may prove this result in a similar way to that of [17, Lemma 9.3 p. 770] where an L^ 
estimate on the translations is proven. Indeed, the proof of Lemma 9.3 [17] holds in the case p = 1 
considered here for a general partition, while for p > 1, it requires the orthogonality condition satisfied 
by the admissible meshes of [17, Definition 9.1 p 762]. We give here a simpler proof dedicated to the 
case p = 1, using the space, as in Lemma 5.1. 



Let u G C^iM."^). For x,y e R'^, one has: 

\u{x + y) — u{x)\ = \ / V u{x + ty) ■ ydt\ < \y\ / \'Vu{x + ty)\dt. 
Jo Jo 

Integrating with respect to x and using Fubini's Theorem gives the well-known result 

„ d 

||u(- -Fy) -nll^imd) < |y| / \Vu\dx <\y\S2\\Diu\\Limd), (79) 

where Vu = {Dm, . . . , Dau). By density of C^{R'^) in W^^'^{R'^), Inequality (79) is also true for 
u G VF^'i(M'^). 

We proceed now as in Lemma 5.1, using the same notations. Let u G BV and Un = u-k p^. Since 
Un G VF^'^(]R"), Inequality (79) gives, for ah y G M'', ||n„(- + y) - 'UnllLi(Rd) < \y\ Ya=i II A'UnllLi(Rd)- 
But, for i = 1, . . . , d, as in Lemma 5.1, H-DjUnH^i^jgd^ < Then, since Un ^ u in L^{R'^), as 

n — > 00, we obtain: 

d 

\\u{- + y)-u\\Li(M.i)<\y\^\\Diu\\M = \y\\\u\\BV yu£ BV, \/y gR'^. (80) 

1=1 

Let u G H^{Q,). One sets u = outside Q so that u G L^(M'^); thanks to lemma 5.1, 11^11^^ ^ 
V^lliilli.i.AI and thus: 

||n(- + y) - u||Li(i;d) < |y| Vd||u||i,i,A4 Vy G R'^. 
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□ 

An easy consequence of Lemmas 5.1 and 5.4 is a compactness result in given in the following 
lemma. 

Lemma 5.5 Let d > 1 and let Q be an open bounded set in M.'^, such that its boundary dO, is a finite 
union of parts of hyperplanes. Let T be a family of polyhedral partitions ofQ in the sense of Definition 
5.1. For M. £ T, let um G H_m{Q) and assume that there exists C G M such that for all A4 G T, 
\Wm\\i,i,M ^ C. Then, the family {uM)MeT is relatively compact in L^{Vt) and also in L^{W^) taking 
Um = outside il.. 

Proof. By Lemma 5.1, the family (n_A/()_A/(gi? is bounded in L^ (Q). Since 17 is bounded, the 
family {u_m)_\4^p is bounded in L^{D.) and also in L^(M'^), taking um = outside $7. Thanks to the 
Kolmogorov compactness theorem, Lemma 5.4 gives that the family {u_M)Me:F is relatively compact 
in L^[Q) and also in L^(M'^) taking um = outside Q. □ 

Note that in fact, the above result also holds for general (non polyhedral) partitions of fi, for instance 
in the case of curved boundaries. In the case p > 1, we need an additional hypothesis on the meshes 
which we state in the following lemma. 

Lemma 5.6 Let d>l, l<p<oo and Q. be a polyhedral open bounded connected subset ofM.'^. Let 
F be a family of meshes of ft in the sense of Definition 2.1. Let rj > be such that, for all D € F, 
one has rj < dK,a/dL,(T < 1/f? for all a & £, where Ai^ = {K,L}. For V G F, let ut> G Hx>{^) and 
assume that there exists C G M such, for all D ^ F , ||ttD||i,p,Al < C . Then, the family {ux))v&F is 
relatively compact in LP[U) and also in LP{M/^) taking u-p = outside Q.. 

Proof. Thanks to Lemma 5.3 and to the fact that $7 is bounded, the family {uxi)v&F is bounded 
in L^{Vl) and also in L^{W^) taking ux) = outside Vl. Thanks once again to the fact that $7 is 
bounded, the family (||tiD||i,i,A^)r>G-F is bounded in M. Then, as in the previous lemma, the Kolmogorov 
compactness theorem gives that the family {ux>)T>eF is relatively compact in L^(f2) and also in L^(M'^) 
taking ut> = outside $7. 

In order to conclude we use, once again. Lemma 5.3. It gives that the family {ut>)v£F is bounded in 
L'^{Q) for some q > p. With the relative compactness in L^{Q,), this leads to the fact that the family 
{ux>)veF is relatively compact in LP{Q) (and then also in L'p{W^) taking uj) = outside J7). □ 



5.2.2 Regularity of the limit 

With the hypotheses of Lemma 5.6, assume that uxi — > u in L'^ as size(P) (Lemma 5.6 gives that 
this is possible, at least for subsequences of sequences of meshes with vanishing size). We prove below 
that u G W^'^in). 

Lemma 5.7 Let d>l, l<p<oo and let Q, be a polyhedral open bounded connected subset ofW^. 
Let {Dn)nGN be a family of discretisations of ft in the sense of Definition 2.1. Let i] > be such 
that, for any discretisation Vn = (A^„, "P^); one has rj < dK,a/dL,a < 1/^ for all a £ £, where 
= {K,L^. For n G N, let li^"^ G Hxi„{i^) and assume that there exists C G M such, for all n G N, 
ll^^"^ ||i,p,Al„ ^ C- Assume also that size('D„) — > as n ^ co. Then: 

1. There exists a sub-sequence of {u^''^^)neN, still denoted by (ii'^"^)„gN, and u G L^(r2) such that 
u*-"-* u in L'P{Q) as n ^ oo. 

2. ue W^'^{n) and 

P-i 

n + ?y)(i p 

l|Vu||ip(Q)d = II |Vu| \\LP{n) < C (81) 

(recall that |Vu| is the Euclidean norm ofVu). 
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Proof. The fact that there exists a subsequence of (ti^"^)nGN) stih denoted by {u^"''')nmi and 
u E LP{fl) such that u^"^ —>■ u in LP{Q) as n ^ oo is a consequence of the relative compactness of 
{u^"''')n£N in LP given in Lemma 5.6. There only remains to prove that u G WQ'P{n). 

Letting u^^^ = and u = outside fi, one also has u^") — > n in L^(M'^). Let us now construct an 
approximate gradient, denoted by bounded in equal to outside Q and converging, 

at least in the distributional sense, to Vu. 

Step 1 Construction of Vpti, for u £ //x)(^)i and its properties. 

Let n G N and V = For this step, one sets u = n^") (not to be confused with the limit of the 
sequence {u^^^)n£f>i)- For a € £, one sets Uo- = if cr is on the boundary of $7. Otherwise, one has 
M-a = {-f^! -f'} and we choose a value between uk and ul (it is possible to choose, for instance, 
Ua = ^{uk + Ul) but any other choice between uk and ul is possible). Then, one defines Vpu on 
K £ D in the following way: 

V-T-ll/. = 

\K 



Vdu = ^ \a\nK,a{ua - uk)- 



The function Vpu is constant on each K £ M and, on using Holder's inequality 
Since X^^-g^;^^, |o"|'i_R',(7 = d\K\, one deduces 



This gives an LP- estimate on Vvu in {LP{Q))'^ (or in {LP{W^))'^, setting Vi^ti = outside Q), in terms 
of ||^i||i,p,Ai, namely 

p-i 

\\\'^vu\\\lv < ^^^"^^"^ " \\u\\i^p,M- (82) 

In order to prove, in the next step, the convergence of this approximate gradient, we now compute 
the integral of this gradient against a test function. Let if G C^(]R'^; M*^), ipK the mean value of (p on 
K £ D, and ip^ the mean value of ip on a. Then, 

/ \7T)U-(fdx= ^ ^ \a\nK,aiua - uk)<Pk = ^ ^ \a\nK,a{-UK)'Pa + R{u,p), (83) 

with 

R{u,(p)= ^ ^ \cr\nK,a{Ua - Uk){PK - ^a)- 

Then, there exists only depending on p, d, p, Q and rj such that \R{u,(p)\ < C;psize(P)||ii||i^p^x. 
Equation (83) can also be written as 

/ Vvu ■ ipdx = / {—Uk) div{ip) dx + R{u,Lp) = — u div{ip) dx + R{u,ip). (84) 

jR'i pr^jy Jk JK.'i 

Step 2 Convergence of Vd^u*^") to Vu and proof of u G Wq'^{Q) . 
We consider now the sequence (n("))„gN- Inequality (82) gives 



IY7 Will ^ i'^ + V)d ^ ||„.(n)|| 
\VVU^ '\\\lp < 11^^^ \\l,p,M- 
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Then, the sequence 

(Vi5n("))„gN is bounded in LP{R'^)'^ and we can assume, up to a subsequence, that 
Vdu^"^ converges to some w weakly in L^(M'^)°', as n — > oo and || \w\ \\lp < I1+!ZM.Z_(7^ 
Let G Cc°°(M'^;M'^), Equation (84) gives 

/ Vpn^") ■ipdx = - [ n(") div(v9) da; + i?(n("\ if). (85) 

Thanks to \R{u^'^\(p)\ < C^size(P„)||n("')||i,p,x„, one has R{u^"'\ip) 0, as n ^ oo. Since u^") ^ u 
in L^(M'^) as n — > oo, passing to the limit in (85) gives: 



/ w ■ (fdx = — u div((^) da;. 



Since ip is arbitrary in C^(M , M ), one deduces that Vu = w. Then u £ W'^'P{R'^) and |||Vti|||LP < 
p-i 

(i+?7)d p ^ Finally, since u = outside Q, one has u £ Wq'^{^1). □ 

6 Conclusion and perspectives 

A symmetric discretisation scheme was introduced for anisotropic heterogeneous problems on dis- 
torted nonconforming meshes. Although this scheme stems from the finite volume analysis, which was 
developed these past years, its formulation is actually derived from a discrete weak formulation; in 
this respect it may be seen as a nonconforming finite element method. Tools of functional analysis 
were obtained, which allow a mathematical analysis of the scheme; the convergence of the discrete 
solution to the exact solution of the continuous problem is shown with no regularity assumption on 
the solution (other than the natural assumption that it is in Hq{^1)). Even though this convergence 
result yields no rate of convergence, it is probably more interesting than error estimates which require 
some assumptions on the diffusion tensor. Nevertheless, we show an order 1 estimate in the case of 
the Laplace operator, which is readily adaptable to regular (say piece- wise C^) isotropic diffusion op- 
erators. The numerical results presented here show the good performance of the scheme (in particular 
order 2 is obtained for the convergence in the norm of the solution), and so do three dimensional 
experiments which were performed in [12] for the incompressible Navier-Stokes equations on general 
grids. Note that the convergence analysis which is performed here readily extends to the non-linear 
setting of Leray-Lions operators. This will be the subject of a future paper. 
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